Uncorrelated Measurements of the Cosmic Expansion History and Dark Energy from 

Supernovae 
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We present a method for measuring the cosmic expansion history H(z) in uncorrelated redshift 
bins, and apply it to current and simulated type la supernova data assuming spatial flatness. If 
the matter density parameter Q m can be accurately measured from other data, then the dark 
energy density history X(z) = px(z) / px(0) can trivially be derived from this expansion history 
H(z). In contrast to customary "black box" parameter fitting, our method is transparent and easy 
to interpret: the measurement of H(z)~ x in a redshift bin is simply a linear combination of the 
measured comoving distances for supernovae in that bin, making it obvious how systematic errors 
propagate from input to output. 

We find the Riess et al. (2004) "gold" sample to be consistent with the "vanilla" concordance 
model where the dark energy is a cosmological constant. We compare two mission concepts for 
the NASA/DOE Joint Dark Energy Mission (JDEM), the Joint Efficient Dark-energy Investigation 
(JEDI), and the Supernova Accelaration Probe (SNAP), using simulated data including the effect of 
weak lensing (based on numerical simulations) and a systematic bias from K-corrections. Estimating 
H(z) in seven uncorrelated redshift bins, we find that both provide dramatic improvements over 
current data: JEDI can measure H(z) to about 10% accuracy and SNAP to 30-40% accuracy. 
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I. INTRODUCTION 

Observational data on type la supernovae (SNe la) in- 
dicate that the expansion of our universe is accelerating 
[1, 2]. This can be explained by the presence of dark en- 
ergy. Various dark energy models have been considered, 
e.g. scalar fields [3-8] and modified gravity [9-15] — see 
[16] and [17] for recent reviews. 

Although current data seem consistent with a cosmo- 
logical constant (e.g., [18-26]), the uncertainties are large 
and more exotic models are not ruled out (e.g., [27- 
31]). To uncover the nature of dark energy, and differ- 
entiate among various dark energy models, it is impor- 
tant that we extract dark energy constraints in a model- 
independent manner [32-34]. The perils of model as- 
sumptions and simplified parametrization of dark energy 
have been shown in [25, 35, 36]. 

Throughout this paper, we assume spatial flatness as 
motivated by inflation. Calibrated cosmological standard 
candles such as SNe la measure the luminosity distance 
cIl(z) = (1 + z)r(z), where the comoving distance 
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(1) 



/0 H(z') * J h(z') 
where c= 1, if* = 100km s _1 Mpc _1 and 

h(z) = H(z)/H* = h(Q) [n m (l + zf + n x X(z)] 

(2) 

with X(z) = px(z)/px(0) denoting the dark energy den- 
sity function. Determining if and (if so) how the dark 
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energy density X(z) depends on cosmic time is the main 
observational goal in the current quest to illuminate the 
nature of dark energy. Given a precise measurement of 
the matter density fraction Q m (from galaxy redshift sur- 
veys, for example), the dark energy density function X(z) 
can be trivially determined from H(z) via equation (2). 

Numerous methods for this have been developed and 
applied in the recent literature, either by parametriz- 
ing X(z) in terms of an equation of state and perhaps 
additional parameters or by aiming for more model- 
independent constraints (e.g., [25]). Equation (1) shows 
that the data directly constrain the cosmic expansion his- 
tory H(z) — in principle. The information-theoretically 
minimal error bars on H (z) that can be obtained from 
supernovae were derived in [37] using a Fisher matrix 
approach. Yet no optimal method for doing this in prac- 
tice has been found other than the "black box" approach 
of parametrizing H (z) somehow and fitting to the data. 
Ideally one would like to measure H(z) in many red- 
shift bins with uncorrelated error bars, but [37] found 
that parametrized fits tend to yield broad and difficult- 
to-interpret window functions, i.e., the measurement in 
a given redshift bin depended also on supernova data far 
outside that redshift range. [29] strengthened this conclu- 
sion by showing that uncorrelated measurements of the 
expansion history (computed by diagonalizing the Fisher 
matrix) tended to probe a broad redshift range. The ap- 
proach most similar to ours is that of [20, 21, 37], where 
H(z)^ 1 is measured by numerically differentiating equa- 
tion (1). These papers tackle the challenge of differenti- 
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FIG. 1: Data used. The measured dimcnsionless comoving dis- 
tances are shown for the Riess et al. (2004) "gold" sample before 
(top panel) and after (bottom panel) flux averaging. The cosmic 
expansion h(z) is simply the inverse of the derivative of this, so the 
visually obvious fact that the derivative is falling immediately tells 
us that the universe was denser and faster expanding in the past. 



ating sparse and noisy data by performing a polynomial 
fit for r(z): [37] performs a global fit whereas [20, 21] 
improve this with a local fit within a sliding window. 

The purpose of the present paper is to solve this prob- 
lem optimally, presenting a method giving uncorrelated 
measurements of the expansion history in arbitrary red- 
shift bins. We will see that this method is both easy 
to implement and easy to interpret. We describe our 
method in Sec. 2. We present our results in Sec. 3 and 
discuss our conclusions in Sec. 4. 



Let us write the comoving distance measured from the 



i th SN la as 



n = r(zi) + rii, 



(4) 



where the noise vector satisfies (n,-) = 0, (n^fij) = crfSi 



A. Transforming to measurements of H(z) 1 

As a first step in our method, we sort the supernovae by 
increasing redshift z\ < Z2 < and define the quantities 
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ft + (n i+ i - rn)/Azi, 



(5) 



where Az, = Zi + \ — z, and f i is the average of 1/H(z) 
over the redshift range (zj, Zi+i). Note that Xi gives an 
unbiased estimate of the average of 1/H(z) in the redshift 
bin, since (x-i) = f t , so the quantities direct (but 

noisy) probes of the cosmic expansion history. Assem- 
bling the numbers Xi into a vector x, its covariance ma- 
trix N = (xx*) — (x)(x)* is tridiagonal, satisfying Nij = 
except for the following cases: 
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The new data vector x clearly retains all the cosmologi- 
cal information from the original data set (the comoving 
distance measurements r,), since the latter can trivially 
be recovered from x up to an overall constant offset by in- 
verting equation (5). In summary, the transformed data 
vector x expresses the SN la information as a large num- 
ber of unbiased but noisy measurements of the cosmic 
expansion history in very fine redshift bins, correspond- 
ing to the redshift separations between neighboring su- 
pernovae. 



B. Averaging in redshift bins 



II. METHOD 

Assuming that the redshifts of SNe la are accurately 
measured, we can neglect redshift uncertainties, and sim- 
ply treat the measured comoving distances r(zi) to the 
supernovae (plotted in Figure 1) as the observables. In 
terms of fiQ, the distance modulus of SNe la, we have 
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1 Mpc 2997.9(1 + z) 
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The second step in our method is to average these noisy 
measurements x into minimum-variance measurements 
t/i, of the expansion history in some given redshift bins. 
For instance, the middle panel of Figure 2 shows an ex- 
ample with seven bins, b = 1,...,7. Let the vector x & 
denote the piece of the x-vector corresponding to the 6 th 
bin, and let denote the corresponding covariance ma- 
trix. Our weighted average yb can then be written 
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FIG. 2: The cosmic expansion history (the dimensionless Hubble 
parameter h(z)) is measured in uncorrelated redshift bins from the 
Riess et al. (2004) "gold" sample (top panel) and from simulated 
future data (middle panel) for the NASA/JDEM mission concepts 
JEDI (solid points) and SNAP (dotted points). The measured 
h(z)~ 1 in a given redshift bin is simply the sum of the comoving 
supernova distances in that bin, weighted by the corresponding 
solid curve in the bottom panel, which roughly speaking subtracts 
more nearby supernovae from more distant ones. 



for some weight vector w b whose components sum to 
unity, i.e., w\ — \ or equivalently e • w fc = 1, where e 
is a vector containing all ones; e< = 1. To find the best 
weight vector w b , we minimize the variance 



Ay b 2 ^(y b 2 )-( 2/b ) 2 =w b 'N b w b 



(8) 



subject to the constraint that the weights add up to unity, 
i.e., that e • w b = 1. This constrained minimization 
problem is readily solved with the Lagrange multiplier 
method, giving 



(9) 



and substituting this back into equation (8) gives the size 
of the corresponding error bar: 



Ay b = (e'N^e) 
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The bottom panel in Figure 2 shows the seven weight 
vectors w & (dotted) corresponding to the seven measure- 
ments yb in the middle panel (solid points). We will also 
refer to the weight vectors as window junctions, since 
they show the contributions to our measurements from 
different rcdshifts. Note that each window function van- 
ishes outside its redshift bin, and that all seven of them 
share a characteristic bump shape roughly corresponding 
to an upside-down parabola vanishing at the bin end- 
points. To illustrate the z-range that each measurement 
probes, we plot it at the median of the window function 
with horizontal bars ranging from the 20th to the 80th 
percentile. 

Since the measurements yb are linear combinations of 
the Xi which are in turn linear combinations of the r%, 
we can also reexpress our measurements directly as lin- 
ear combinations of the original supernova comoving dis- 
tances: 
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would be essentially the negative derivative of the old 
window functions if all redshift intervals were the same. 
These new window functions are also plotted in the bot- 
tom panel of Figure 2 (solid "sawtooth" curves), and are 
seen to be roughly linear (as expected for a parabola 
derivative) , effectively subtracting supernovae at the near 
end of the bin from those at the far end. 

A major advantage of this method is its transparency 
and simplicity. If one fits some parametrized model of 
H(z) to the SN la data by maximizing a likelihood func- 
tion, then the resulting parameter estimates will be some 
complicated (and gencrically nonlinear) functions of all 
the data points fj. In contrast, the measurement y§ in 
our 5 th bin in figure 2 (middle panel) is simply a linear 
combination of the comoving distance measurements for 
the supernovae in the 5 th bin (1.0 < z% < 1.4) as defined 
by the 5 th window function in the bottom panel, so it 
is completely clear how each particular supernova affects 
the final result. In particular, the supernovae outside of 
this redshift range do not affect the measurement at all. 

We conclude this section by discussing some details 
useful for the reader interested in applying our method 
in practice. 



C. Creating uncorrelated redshift bins 

We suggest discarding those Xi straddling neighboring 
bins, i.e., whose two supernovae fall on either side of a bin 
boundary. For say 7 bins there are only 6 such numbers, 
so this involves a rather negligible loss. The advantage 
is that it ensures that two measurements yb and y^ have 
completely uncorrelated error bars if b ^= b' , since their 
window functions w have no supernovae in common. 
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D. Flux averaging 

To minimize the bias in the iJ (z)-measurement due to 
weak lensing, we use flux-averaging [38, 39] 1 . Specifically, 
we compress the full supernova data set T\ into a smaller 
number of of flux-averaged supernovae assigned to the 
mean redshift in each of a large number of bins of width 
Az. We use Az = 0.05 for the current data and Az = 
0.005 for the simulated data. Note that since we have 
assumed a Gaussian distribution in the magnitudes of 
SNe la at peak brightness, flux-averaging leads to a tiny 
bias of —crf nt In 10/5 mag [41]. We have removed this 
tiny bias in the data analysis. 

As a side effect, this averaging in narrow bins makes 
the denominators Azj roughly equal in equation (5), so 
that the only noticeable source of wiggles in the window 
functions in Figure 2 (bottom) is Poisson noise, i.e., that 
some of these narrow bins contain more supernovae than 
others. The method of course works without this aver- 
aging step as well. In that case, the window functions w 
wiggle substantially because of variations in the redshift 
spacing between supernovae, since very little weight is 
given to Xi if Azi happens to be tiny. However, we find 
that the supernova window functions w remain rather 
smooth and well-behaved functions of redshift, as ex- 
pected — two supernovae very close together with the 
same noise level Oi automatically get the same weight. 
This means that our method effectively averages such 
similar redshift supernovae anyway, even if we do not do 
so by hand ahead of time. The difference between flux 
averaging and this automatic averaging is simply that we 
average their fluxes rather than their comoving distances 
— these two types of averaging are not equivalent since 
the flux is a nonlinear function of the comoving distance. 



III. RESULTS 

Figure 2 shows the results of applying our method to 
both real data (top panel) and simulated data (middle 
panel) , with the dimensionless expansion rate of the Uni- 
verse h(z) measured in between three and seven uncor- 
related redshift bins. 

The top panel uses the "gold" set of 157 SNe la pub- 
lished by [40]. The error bars are seen to be rather large, 
and consistent with a simple flat Q m = 0.3 concordance 
model where the dark energy is a cosmological constant. 2 
As was shown in [37] using information theory, the rel- 
ative error bars on the cosmic expansion history H(z) 
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FIG. 3: Same as Fig. 2(b), but with a systematic bias of dm aya = 
0.02 z from K-correction uncertainties added in addition to lensing 
noise computed from numerical simulations (Barber 2003 [42]). 



scale as 



H iW2A z 3/2' ^"J 

for N supernovae with noise o. Here Az is the width 
of the redshift bins used, so one pays a great price for 
narrower bins: halving the bin size requires eight times 
as many supernovae. The origin of this (Az)~ 3 / 2 -scaling 
is intuitively clear: the noise averages down as (Az) -1 / 2 , 
and there is an additional factor of (Az) -1 from effec- 
tively taking the derivative of the data to recover H(z)~ 1 
from the integral in equation (5) 3 . The bottom panel in 
Figure 2 shows that the method effectively estimates this 
derivative by subtracting supernovae at the near end of 
the bin from those at the far end of the bin and dividing 
by Az. 

This means that for accurately measuring H(z) and 
thereby the density history of the dark energy, numbers 
really do matter. For example, in order to measure the 
dark energy density function to 10% accuracy in seven 
uncorrelated redshift bins as in the middle panel of Figure 
2 (with Az = 0.243), we need to have around 14,000 SNe 
la. We have simulated SN la data by placing supernovae 
at random redshifts, with the number of SNe la per 0.1 
redshift interval given by a distribution. The intrinsic 
brightness of each SN la at peak brightness is drawn from 
a Gaussian distribution with a dispersion of (Ji n t = 0.16 
mag. As the fiducial cosmological model, we used a flat 
universe with J7a = 0.7. 

We compare two mission concepts for the NASA/DOE 
Joint Dark Energy Mission (JDEM): the Joint Efficient 
Dark-energy Investigation ( JEDI) [43] and the Supernova 
Accelaration Probe (SNAP) [44]. For JEDI (solid points 
in Figure 2), the number of SNe la per 0.1 redshift in- 
terval is obtained by fitting the measured SN la rate as 



1 A Fortran code that uses flux-averaging statistics to compute 
the likelihood of an arbitrary dark energy model (given the SN 
la data from [40]) can be found at http : / /www. nhn.ou.edu/ ~ 
wang / S N code / . 

2 This is consistent with the findings of [25]. 



Anaogous estimates of the equation of state w(z) have a painful 
(Az) - 5 / 2 -scaling, since they effectively involve taking the second 
derivative of the data. 
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function of redshift [45-47] to a model assuming a con- 
servative delay time between star formation and SN la 
explosion of 3.5 Gyrs. For SNAP (dotted points in Fig- 
ure 2), the number of SNe la per 0.1 redshift interval is 
taken from Figure 9 in [44]. 

We consider two kinds of SN la systematic uncer- 
tainties: weak lensing due to intervening matter and 
a systematic bias due to K-corrections. We include 
the weak lensing effect by assigning a magnification /i 
drawn from a probability distribution p(p), extracted 
using an improved version of the Universal Probability 
Distribution Function (UPDF) method [48] from the nu- 
merical simulations of weak lensing by Barber (2003) 
[42]. The total uncertainty in each SN la data point is 
V a int + a iens( z ) 2 , with <Ji ens (z) extracted from Barber 
(2003) [42]: 

ai ens (z) ~ 0.00311 + 0.08687z - 0.00950z 2 (14) 

We consider a systematic bias of Am sys = 0.02z due to 
K-corrections following [32] . 

We did not include the systematic bias due to K- 
corrections in Fig. 2, in order to compare the real data 
[40] and simulated data on an equal footing. 

In Fig. 3, we show the effect of adding the systematic 
bias due to K-corrections in addition to the weak lensing 
noise. Comparing Fig. 3 with Fig. 2 (b), we see that 
the systematic bias does not have a significant effect on 
the uncorrelated estimates of H(z). This is because our 
method effectively reduces a global systematic bias into 
a local bias with a much smaller amplitude (see Sec. 2). 

Figure 4 shows how the recovery of the cosmic expan- 
sion history h(z) depends on the number of redshift bins 
(6, 7 and 8), assuming no systematic bias and no lens- 
ing, and agrees well with the theoretical (Az) 3 / 2 scaling. 
Contrasting current and future data with roughly the 
same redshift binsize for both (nbj„=5 for current data, 
and nfci„=7 for the simulated data) shows that JEDI 
shrinks the error bars by more than an order of mag- 
nitude, so the potential improvement with a successful 
JDEM would be dramatic. 

Figure 4 compares three different data sets: (a) Half of 
the JEDI data, with a reduced intrinsic scatter of (j\ n t = 
0.08 mag from sub-typing, (b) All the JEDI data, with 
a int = 0.16 mag. (c) All the SNAP data (plus 25 SNe la 
at z < 0.1), with ai n t = 0.16 mag. 

Note that since the measured quantity h(z)^ 1 typically 
is not a straight line, the measured average of this curve 
over a redshift bin will generally lie either slightly above 
of below the curve at the bin center. Figure 2 shows that 
this bias is substantially smaller than the measurement 
uncertainties, since h(z) and /i(z) _1 are rather well ap- 
proximated by straight lines over the narrow redshift bins 
that we have used. 

A second caveat when interpreting our figures is that 
the absolute calibration of SN la is not perfectly known 
— changing this simply corresponds to multiplying the 
function h(z) by a constant, i.e., to scaling the measured 
curve vertically. 
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FIG. 4: How the recovery of the cosmic expansion history h(z) 
depends on the number of redshift bins, assuming no systematic 
bias and no lensing. (a) Half of JEDI data, with a reduced intrinsic 
scatter of <r; nt = 0.08 mag. (b) All of JEDI data, with cri nt = 0.16 
mag. (c) All of SNAP data (plus 25 SNe la at z < 0.1), with 
o-int = 0.16 mag. 

IV. DISCUSSION 

We have presented a method for measuring the cosmic 
expansion history H(z) in uncorrelated redshift bins, and 
applied it to current and simulated supernova la data 
assuming spatial flatness. Whereas previously proposed 
approaches involve "black box" parameter fitting, this 
method is transparent and simple to interpret: the mea- 
surement of H(z)^ 1 in a redshift bin is simply a linear 
combination of the measured comoving distances for su- 
pernovae in that bin, with weights that roughly corre- 
spond to subtracting closer supernovae from more dis- 
tant ones. Such transparency is particularly helpful for 
understanding how systematic errors in the input affect 
the output. For instance, a constant systematic dimming 
throughout a redshift bin will leave that measurement 
unaffected. 

This method is useful for understanding the nature of 
dark energy, since the dark energy density history follows 
from this expansion history H(z) if the matter density 
parameter Q m can be accurately measured from other 
data (e.g., the cosmic microwave background, galaxy 
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clustering and gravitational lensing). We found the Ricss 
et al. (2004) "gold" sample to be consistent with the 
"vanilla" concordance model where the dark energy is 
a cosmological constant, but that much larger numbers 
of supernovae are needed for true precision tests of the 
nature of dark energy. We obtain good agreement with 
the results of [20, 21] — the fact that two quite different 
techniques give consistent answers indicates a reassuring 
robustness to method and data details. 

Looking towards the future, we compare two mission 
concepts for the NASA/DOE Joint Dark Energy Mission 
(JDEM), the Joint Efficient Dark-energy Investigation 
(JEDI) and the Supernova Accelaration Probe (SNAP), 
using simulated data including the effect of weak lens- 
ing and bias from K-corrections. Estimating H{z) in 
seven uncorrelated redshift bins, we find that both pro- 
vide dramatic improvements over the present state-of-the 
art: JEDI can measure H(z) to about 10% accuracy, and 
SNAP can measure H{z) to 30-40% accuracy (Figure 2). 

Our results show that numbers do matter. For exam- 
ple, in order to measure the cosmic expansion history to 
10% accuracy in seven uncorrelated redshift bins (with 
Az = 0.243), we need to have around 14,000 SNe la 
(as expected from two years of JEDI data), assuming a 
dispersion of 0.16 magnitudes in SN la peak brightness 
(Figure 2). Estimating H(z) to 10% in smaller uncor- 
related redshift bins will require an even larger number 
of SNe la (if a significant reduction in intrinsic disper- 
sion is not assumed for a sizable fraction of the SNe la), 



which will be difficult to obtain in a feasible two-year 
space mission. Also, H(z) estimated in Az = 0.243 red- 
shift bins (for < z < 1.7) to 10% accuracy will give 
us a powerful means to differentiate between a cosmo- 
logical constant and dark energy models which are not 
fine-tuned to mimic a cosmological constant. A sample 
with a large number of SNe la allows tighter calibration 
of SNe la as standard candles and subtyping of SNe la to 
reduce diversity. This may yield a smaller set of SNe la 
with substantially smaller intrinsic dispersion, which can 
lead to more robust and stable estimates of the expan- 
sion history of the universe (Figure 4a). If the subtyping 
works well in reducing the intrinsic dispersion of SNe la, 
we can expect to be able to measure H{z) in more than 
seven redshift bins to 10% accuracy for < z < 1.7 (Fig- 
ure 4a). 

A precise measurement of the cosmic expansion his- 
tory as a free function of cosmic time to 10% accuracy 
would represent a dramatic improvement in our knowl- 
edge about dark energy. Our results suggest that a 
JDEM can achieve this scientific goal. 
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